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Abstract 

We present the results of a classical molecular dynamic simulation as well as of an ab initio molecular dynamic 
simulation of an amorphous silica surface. In the case of the classical simulation we use the potential proposed by 
van Beest et al. (BKS) whereas the ab initio simulation is done with a Car-Parrinello method (CPMD). We find 
that the surfaces generated by BKS have a higher concentration of defects (e.g. concentration of two-membered 
rings) than those generated with CPMD. In addition also the distribution functions of the angles and of the 
distances are different for the short rings. Hence we conclude that whereas the BKS potential is able to reproduce 
correctly the surface on the length scale beyond «5 A, it is necessary to use an ab initio method to predict reliably 
the structure at small scales. 
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1. Introduction 

Obtaining a good understanding of the struc- 
tural and dynamical properties of the surface of 
amorphous silica is very important for the manu- 
facture of glass as well as the construction of elec- 
tronic devices [1] . This is the reason why in the past 
a large number of experiments have been done to 
investigate this type of surface. Since in real exper- 
iments it is rather difficult to obtain reliably details 
on the structure also quite a few computer simu- 
lations have been done in order to study this sys- 
tem (see [2] and references therein). Most of these 
studies have, however, been done by using effec- 
tive classical potentials, such as, e.g., the one pro- 



posed some years ago by van Beest, Kramer, and 
van Santen (BKS) [3]. Although it has been shown 
that these type of potentials can reproduce quite 
reliably the structure and dynamics of silica in the 
bulk ( see, e.g. , [4] and references therein) , it is much 
less obvious to what extent they are also able to 
give a correct description of the properties of silica 
close to a surface, since the parameters for these 
potentials, effective charges, etc., have often been 
optimized to reproduce only experimental data for 
the bulk. One possibility to avoid this problem with 
the classical effective potential is to use ab initio 
simulations such as the scheme proposed by Car 
and Parrinello [5] since in this type of approach an 
effective potential between the ions is calculated 
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self consistently on the fly, i.e. the instantaneous 
geometry of the ions is always taken into account. 
The drawback of this approach is that due to the 
huge computational burden only relatively short 
time scales, a few ps, as well as small systems, a 
few hundred particles, can be simulated, whereas 
classical simulations allow to simulate thousands 
of particles over several ns. 

In the present work we compare the results of a 
classical simulation of an amorphous silica surface 
with the BKS potential with the results obtained 
by the Car-Parrinello Molecular Dynamics method 
(CPMD). The goal is to check which quantities are 
reproduced correctly by the BKS potential, using 
the results of the CPMD simulation as the refer- 
ence system. 



2. The BKS potential and the setup of the 
geometry 

We have first prepared the system using the BKS 
potential. In this two-body potential the atoms in- 
teract also by means of a Coulomb potential where 
the effective charges of a silicon and oxygen atom 
is 2.4 and —1.2, respectively. More details on this 
potential can be found in Ref. [3] . 

In order to minimize finite size effects as well 
as surface effects it is customary to use periodic 
boundary conditions (PBC) in all three directions. 
If one wants to investigate a free surface, the most 
straightforward idea is to use a fllm geometry, i.e. 
to have PBC in two directions and to have an infi- 
nite free space above and below the system. Unfor- 
tunately it turns out, however, that from a compu- 
tational point of view this setup is not very good 
for systems with Coulombic interactions (i.e. such 
as the one studied here) , since it prevents to make 
an efficient use of the Ewald summation method. 
Therefore we have adopted the following strategy 
which is explained also in Fig. 1 (see also [6,7]): i) 
We start with a relatively large (bulk) silica system 
at T = 3400 K with PBC in three dimension (with 
box size = Ly = 11. 51 A and L'^ = 23A. ii) 
We cut the system perpendicular to the z-direction 
into two pieces. Without loss of generality we can 
assume that the mean position of this cut is at 
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Fig. 1. Procedure to create the used sandwich geometry 

z — L'^. At this point it is very important, that 
we cut only oxygen-silicon-bonds, such that we get 
only free oxygen atoms at this interface (thus the 
interface will have a bit of roughness), iii) These 
free oxygen atoms are now saturated by hydrogen 
atoms. The place of these hydrogen atoms are cho- 
sen such that each of the new oxygen-hydrogen 
bonds is in the same direction as the oxygen-silicon 
bond which was cut and has a length of approx- 
imately 1 A. The interaction between the hydro- 
gen atoms and the oxygen atoms as well as the 
silicon atoms are described only by a Coloumbic 
term. The value of the effective charge of the hy- 
drogen atoms is set to 0.6, which ensures that the 
system is still (charge) neutral, iv) We make atoms 
which have a distance from this interface that is 
less than 4.5 A completely immobile, whereas the 
atoms that have a larger distance can propagate 
subject to the force field, v) We add in z-direction 
an empty space of Az = 6.0 A and thus generate a 
free surface at around 14.5 A. With this sandwich 
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geometry we now can use periodic boundary con- 
ditions in all three directions. We have made sure 
that the value of Az is sufficiently large that the 
results do not depend on it anymore [8] . Note that 
it is not advisable to choose Az too large, since this 
would increase the cost of the CPMD simulation. 
At the end of this procedure we have a system of 91 
oxygen, 43 silicon and 10 hydrogen atoms in a sim- 
ulation box with La;=L„=11.51 A and « 25 A. 



3. CPMD-Simulation 

Since the time scale which is accessible to the 
CPMD-mcthod is very restricted, we have to 
combine the ab initio with classical calculations. 
For this we first prepared a classical system as 
described in the previous section, equilibrated it 
for about 1 ns which is sufhcient to equilibrate 
it completely [4]. From a subsequent production 
run with the same duration we picked 100 statis- 
tical independent configurations and used them 
to characterize the static properties of the sys- 
tem with a high accuracy. Using a subset of these 
configurations as starting points, we subsequently 
started CPMD-simulation using the CPMD code 
developed in Stuttgart [9]. For the CPMD we 
used conventional pseudopotentials for silicon and 
oxygen and the BLYP exchange-functions [10,11]. 
The electronic wave-functions were expanded in a 
plane wave basis set with an energy cutofi^ of 60 Ry 
and the equations of motion were integrated with 
a time step of 0.085 fs for 0.2 ps. In the analysis 
of the CMPD data only those configurations were 
taken into account that were produced later than 
5 fs aftca- the start of the CPMD run in order to 
allow the system to equilibrate at least locally [12]. 

In the analysis of the classical configurations we 
noted that typically one of the three following sit- 
uations is present on the surface: 

- systems with no defects (i.e. all Si and O atoms 
are four and two-fold coordinated, respectively) 

- systems with an undercoordinated oxygen atom 
and an undercoordinated silicon atom 

- systems with an overcoordinated oxygen atom 
and an undercoordinated oxygen atom 
Therefore we picked for each case two BKS con- 



figurations and started the CPMD runs. 

The largest differences between the results of the 
classical and of the CPMD simulation are found 
for the short rings (n<5). (A ring is a closed loop of 
n consecutive Si-0 segments [13].) Figure 2 shows 
the probability to find a ring of size n for the case 
of BKS and CPMD. We see that the BKS potential 
overestimates the frequency with which a ring of 
size two occurrs by about a factor of two. Related 
to this is the observation that the overshoot that is 
observed in the z— dependent mass density profile 
is less pronounced in the case of the CPMD than 
the one for the BKS (inset of Fig. 2), since two- 
membered rings are relatively dense. 
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Fig. 2. Probability to find a ring of size n. Inset: 
z— dependence of the mass density. 

Prom this inset we also see that the density pro- 
file for the CPMD extends to larger z values that 
the one for the BKS (by about 0.4 A). This is be- 
cause the BKS potential is not able to reproduce 
correctly the density of silica at zero pressure. 

Another interesting result is the dependence of 
the distribution of angles 0-Si-O on the ring size 
(Fig. 3). For large n, n >4, i.e. the sizes which are 
normally found in the bulk [14], the results of the 
two different methods are in good agreement [12]. 
For smaller n. however, the mean 0-Si-O-angle 
from CPMD is shifted to larger values in compar- 
ison to the classical one. This shift becomes more 
pronounced with decreasing n. Furthermore also 
the shape of the distributions starts to become dif- 
ferent if n is small. 

This effect can be understood better by analyz- 
ing the partial radial distribution functions g{r) 
which are shown in Fig. 4. We see that for the Si-0 
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Fig. 3. Distribution of O-Si-O-angles for different ring sizes. 
BKS=filled symbols; CPMD=open symbols. 

pair the curves from CPMD are shifted to larger 
distances by about 0.04 A and that this shift is in- 
dependent of n. Also the g{r) for the 0-0 pairs are 
shifted to larger r, but this time the amount does 
depend on n. In particular we note that the 0-0 
distance is nearly independent of n for the case of 
CPMD, whereas it increases with n for the case of 
BKS. These effects results in the difference in the 
distribution of the 0-Si-O angles if n is small. 




Fig. 4. Radial distribution function for different ring sizes. 
Left: O-O. Right Si-O. BKS=filled symbols; CPMD=open 
symbols. 



4. Conclusion 

In this work we have investigated some struc- 
tural properties of an amorphous silica surface. In 
particular we have studied how these properties de- 
pend on the simulation method: A classical simula- 



tion with the potential proposed by van Beest et al. 
(BKS) and a Car-Parrinello simulation (CPMD). 
We find that the structure on larger length scales 
are independent of the method used, whereas the 
details of structural elements on short scales (short 
rings, distribution function for angles, etc.) differ. 
Thus this shows that it is probably necessary to 
use ah initio methods if one wants to understand 
these systems at short length scales in a quantita- 
tive way. 
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